# -*- mode: python; coding: utf-8; fill-column: 100 -*-

# simpintarith.py -- Naïve implementation of interval arithmetics
# Copyright (C) 2006 Matti Hänninen

# This program is free software; you can redistribute it and/or modify it under the terms of the GNU
# General Public License as published by the Free Software Foundation; either version 2 pof the
# License, or (at your option) any later version.

# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
# even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.

# You should have received a copy of the GNU General Public License along with this program; if not,
# write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
# 02110-1301, USA.

# $version: 0.0.1$
# $author: Matti Hänninen <matti.hanninen@iki.fi>$
# $created: 2006-Jul-29$
# $last-updated: 2006-Jul-30$

"""Naïve implementation of interval arithmetics."""

__all__ = ["Interval", "cont", "iacos", "iasin", "iatan", "iceil",
           "icos", "iexp", "ifloor", "ilog", "interval", "ipow",
           "irecip", "isin", "isq", "isqrt", "itan", "olap", "olapx"]

import math

def _is_numeric(obj):
    return isinstance(obj, (int, long, float))

def _are_numerics(*objs):
    for obj in objs:
        if not isinstance(obj, (int, long, float)):
            return False
    return True

def _is_interval(obj):
    return isinstance(obj, Interval)

def _are_intervals(*objs):
    for obj in objs:
        if not isinstance(obj, Interval):
            return False
    return True

def _unsup_type_error(op, *args):
    types = ["'%s'" % arg.__class__.__name__ for arg in args]
    last = types[-1]
    init = types[:-1]
    if not init:
        msg = "unsupported operand type for %s: %s" % (op, last)
    else:
        msg = "unsupported operand types for %s: %s" % (op, ', '.join(init) + ' and ' + last)
    return TypeError(msg)

def interval(obj):
    if isinstance(obj, Interval):
        return obj
    elif _is_numeric(obj):
        return Interval(obj)
    else:
        raise TypeError("an interval or non-complex numeric required")
        
class Interval:
    def __init__(self, a, b=None):
        if b is None:
            if _is_numeric(a):
                self.l = self.u = a
            else:
                raise ValueError()
        elif _are_numerics(a, b):
            if a < b:
                self.l = a
                self.u = b
            else:
                self.l = b
                self.u = a
        else:
            raise ValueError()

    def __repr__(self):
        return '<Interval [' + str(self.l) + ',' + str(self.u) + '] with id ' + str(id(self)) +'>'

    def __str__(self):
        if self.l == self.u:
            return '[%s]' % str(self.l)
        else:
            return '[%s,%s]' % (str(self.l), str(self.u))

    def __pos__(self):
        return self

    def __neg__(self):
        return Interval(-self.u, -self.l)

    def __abs__(self):
        if self.l > 0:
            return self
        elif self.u < 0:
            return self.__neg__()
        else:
            raise ArithmeticError("Absolute value not well defined")

    def __add__(self, other):
        if _is_interval(other):
            return Interval(self.l + other.l, self.u + other.u)
        elif _is_numeric(other):
            return Interval(self.l + other, self.u + other)
        else:
            raise _unsup_type_error("+", self, other)

    __radd__ = __add__                  # TODO: Type names get presented in wrong order in error
                                        # message. To fix it we may have to write separate function
                                        # for this one too.

    def __iadd__(self, other):
        if _is_interval(other):
            self.l += other.l
            self.u += other.u
        elif _is_numeric(other):
            self.l += other
            self.u += other
        else:
            raise _unsup_type_error("+=", self, other)

    def __sub__(self, other):
        if _is_interval(other):
            return Interval(self.l - other.u, self.u - other.l)
        elif _is_numeric(other):
            return Interval(self.l - other, self.u - other)
        else:
            raise _unsup_type_error("-", self, other)

    def __rsub__(self, other):
        if _is_interval(other):
            return Interval(other.l - self.u, other.u - self.l)
        elif _is_numeric(other):
            return Interval(other - self.u, other - self.l)
        else:
            raise _unsup_type_error("-", other, self)

    def __isub__(self, other):
        if _is_interval(other):
            self.l -= other.u
            self.u -= other.l
        elif _is_numeric(other):
            self.l -= other
            self.u -= other
        else:
            raise _unsup_type_error("-=", self, other)

    def __mul__(self, other):
        if _is_interval(other):
            a = self.l * other.l
            b = self.l * other.u
            c = self.u * other.l
            d = self.u * other.u
            return Interval(min(a, b, c, d), max(a, b, c, d))
        elif _is_numeric(other):
            if other >= 0:
                return Interval(self.l * other, self.u * other)
            else:
                return Interval(self.u * other, self.l * other)
        else:
            raise _unsup_type_error("*", self, other)

    __rmul__ = __mul__                  # TODO: Type names get presented in wrong order in error
                                        # message. To fix it we may have to write separate function
                                        # for this one too.

    def __imul__(self, other):
        if _is_interval(other):
            a = self.l * other.l
            b = self.l * other.u
            c = self.u * other.l
            d = self.u * other.u
            self.l = min(a, b, c, d)
            self.u = max(a, b, c, d)
        elif _is_numeric(other):
            if other >= 0:
                self.l *= other
                self.u *= other
            else:
                self.l, self.u = self.u * other, self.l * other
        else:
            raise _unsup_type_error("*=", self, other)
                                 
    def _reciprocal(self):
        if self.l > 0 or self.u < 0:
            return Interval((1).__truediv__(self.u), (1).__truediv__(self.l))
        else:
            raise ZeroDivisionError("interval division by interval containing zero")
    
    def __div__(self, other):
        if _is_interval(other):
            return self.__mul__(other._reciprocal())
        elif _is_numeric(other):
            if other > 0:
                return Interval(self.l.__truediv__(other), self.u.__truediv__(other))
            elif other < 0:
                return Interval(self.u.__truediv__(other), self.l.__truediv__(other))
            else:
                raise ZeroDivisionError("interval division by scalar zero")
        else:
            raise _unsup_type_error("/", self, other)
    
    __truediv__ = __div__

    def __rdiv__(self, other):
        if _is_interval(other):
            return other.__mul__(self._reciprocal())
        elif _is_numeric(other):
            raise NotImplementedError()
        else:
            raise _unsup_type_error("/", other, self)

    __rtruediv__ = __rdiv__

    def __idiv__(self, other):
        if _is_interval(other):
            a = self.l * other.l
            b = self.l * other.u
            c = self.u * other.l
            d = self.u * other.u
            self.l = min(a, b, c, d)
            self.u = max(a, b, c, d)
        elif _is_numeric(other):
            if other > 0:
                self.l = self.l.__truediv__(other)
                self.u = self.u.__truediv__(other)
            elif other < 0:
                self.l, self.u = self.u.__truediv__(other), self.l.__truediv__(other)
            else:
                raise ZeroDivisionError("interval division by scalar zero")
        else:
            raise _unsup_type_error("/=", self, other)

    def __cmp__(self, other):
        return NotImplemented

    def __rcmp__(self, other):
        return NotImplemented

    def __eq__(self, other):
        return NotImplemented

    def __lt__(self, other):
        return NotImplemented

    def __le__(self, other):
        return NotImplemented

    def __gt__(self, other):
        return NotImplemented

    def __ge__(self, other):
        return NotImplemented

def irecip(interval):
    """irecip(x)

    Returns the reciprocal of the interval x.  If x contains zero, the
    function raises ZeroDivisionError."""
    return interval._reciprocal()

def cont(x, y):
    """cont(x, y)

    Tests whether the interval x contains y.  If y is an interval then
    it is contained when all of its points are also in x.  If y
    is a numeric then it is contained if it is in x."""
    if _is_interval(x):
        if _is_interval(y):
            return x.l <= y.l and y.u <= x.u
        elif _is_numeric(y):
            return x.l <= y and y <= x.u
        else:
            raise _unsup_type_error("'cont'", x, y)
    else:
        raise _unsup_type_error("'cont'", x, y)

def olap(x, y):
    """olap(x, y)

    Tests whether intervals x and y overlap, i.e., they share at least
    one common (end) point."""
    if _are_intervals(x, y):
        if x.u < y.l or y.u < x.l:
            return False
        else:
            return True
    else:
        raise _unsup_type_error("'olap'", x, y)

def olapx(x, y):
    """olapx(x, y)

    Tests whether intervals x and y overlap extensively, i.e., they
    share more than one single common point."""
    if _are_intervals(x, y):
        if x.u <= y.l or y.u <= x.l:
            return False
        else:
            return True
    else:
        raise _unsup_type_error('olapx', x, y)

_twopi = 2.0 * math.pi

def iceil(x):
    """iceil(x)

    Returns the ceiling of the interval x."""
    raise NotImplementedError()

def ifloor(x):
    """ifloor(x)

    Returns the floor of the interval x."""
    raise NotImplementedError()

def iexp(x):
    """iexp(x)

    Returns the (natural) exponential of the interval x."""
    raise NotImplementedError()

def ipow(x, y):
    """ipow(x, y)

    Returns the interval x raised to the y-th power.  Note that only
    integral powers are allowed when x overlaps with negative real axis.
    When x doesn't overlap with y both fractional and interval powers
    are allowed."""
    raise NotImplementedError()

def ilog(x, y):
    """ilog(x, y)

    Returns the base y logarithm of the interval x.  The interval x must
    not overlap with the interval (-inf , -1].  If y is not given then
    it assumed to be e."""
    raise NotImplementedError()

def isq(x):
    """isq(x)

    Returns the square of the interval x."""
    if _is_interval(x):
        if x.l >= 0:
            return Interval(x.l * x.l, x.u * x.u)
        elif x.u <= 0:
            return Interval(x.u * x.u, x.l * x.l)
        else:
            return Interval(0, max(x.l * x.l, x.u * x.u))
    else:
        raise TypeError("an interval is required")

def isqrt(x):
    """isqrt(x)

    Returns the positive square root of the interval x.  The interval x
    must not overlap with the negative real axis."""
    raise NotImplementedError()

def isin(x):
    """isin(x)

    Returns the sine of the interval x."""
    raise NotImplementedError()

def iasin(x):
    """iasin(x)

    Returns the arc sine of the interval x in radians."""
    raise NotImplementedError()

def icos(x):
    """icos(x)

    Returns the cosine of the interval x."""
    if not _is_interval(x):
        raise TypeError("an interval is required")
    l = x.l % _twopi
    u = x.u % _twopi
    cl = math.cos(l)
    cu = math.cos(u)
    if cu < cl: cl, cu = cu, cl
    if l <= u:
        if u <= math.pi or math.pi <= l:
            return Interval(cl, cu)
        else:
            return Interval(-1.0, cu)
    elif l <= math.pi or math.pi <= u:
        return Interval(-1.0, 1.0)
    else:
        return Interval(cl, 1.0)

def iacos(x):
    """iacos(x)

    Returns the arc cosine of the interval x in radians."""
    raise NotImplementedError()

def itan(x):
    """itan(x)

    Returns the tangent of the interval x."""
    raise NotImplementedError()

def iatan(x):
    """iatan(x)

    Returns the arc tangent of the interval x in radians."""
    raise NotImplementedError()

# eof
